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Abstract 

'■^ ' The 2PI effective action formalism for quantum fields out of equilibrium 

is set up in an expanding (Friedmann-Robertson- Walker) background. We 
write down and solve the evolution equations for a ip^ model at O(A^) in 
OO I a coupling expansion. We comment on issues of renormalization, lattice 

*^ ' discretization and the range of applicability of the approach. A num- 

ber of example calculations are presented, including thermalization and 
(p)reheating. Generalizations to more complicated systems and applica- 



■ tions are discussed. 
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1 Introduction 



^ I In recent years, significant attention has been drawn to the process of thermal- 

ization of quantum fields. Quantitative description of the physics of the very 
early Universe and of heavy-ion collision experiments requires an understanding 
of the real-time dynamics of quantum fields at finite energy density, but out of 
equilibrium. 

One very promising development is the application of the 2PI-formalism [H [2l 
|3] , which allows the derivation and explicit numerical solution of a set of equations 
of motion for the mean field and propagator in the full quantum theory. This 
is realized through the truncation of a controlled diagram expansion in terms 
of 2PI diagrams. Already at next-to-leading order (NLO) in either a coupling 
or expansion, interacting systems exhibit equilibration, effective dissipation 
and thermalisation to the quantum equilibrium state [21 HI [5l [6l [TJ [H [9] . 
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The physics of the early Universe is described not in a static background, 
but in expanding space. This is often approximated by a homogeneous, flat 
Friedmann-Robertson- Walker (FRW) space-time, parametrized by the scale fac- 
tor a(t) of the metric ds'^ = dt^ — a?{t)d:K^. In some cases, expansion can be 
neglected when the time scale of a phenomenon is very short on the time scale 
of expansion. But in general, and of course in principle, expansion should be 
included in the description of early Universe physics. 

Many processes at high temperature or energy density are well described by 
the classical approximation, where a Monte-Carlo sample of initial field config- 
urations are evolved using Hamiltonian equations of motion (see [10] for a brief 
review). The observables of interest are then averages over this classical ensem- 
ble. The approximation must however break down eventually, as classical fields 
equilibrate to a classical equilibrium, which suffers from the Rayleigh- Jeans prob- 
lem: in the continuum temperature will go to zero, and on the lattice it will be 
cut-off dependent. 

In this paper, we set out the 2PI formalism in a FRW space-time, and solve 
the resulting equations numerically for some example applications. A number 
of studies have been carried out in this context in the Hartree approximation 
[m [12], [131 [HI HHj (which is also leading order (LO) in a 2PI coupling expansion), 
and even before that, the formalism was set out in [T6lll7j. Recently, attempts 
have been made to partly include the effect of an expanding background for 
specific applications also at NLO [T8l [T9| [20]. 

Two main issues present themselves. Firstly, as we are discretizing the system 
on a finite co-moving lattice, there is only a finite number of momentum modes 
available, and as the lattice expands in time these will be redshifted towards 
the IR in physical units. This means that the physical cut-off changes in time. 
Therefore, there is a limit on how many e-folds one can run the simulation before 
running out of "dynamical range". In practice, this means that at some point 
discretization errors become important, and results can no longer be relied upon. 
As a result, reliably simulating cosmological inflation proper is a daunting task, 
as the Universe expands many e-folds. Still, most of the inflationary stage is often 
well described by semi-analytical tools and the slow-roll approximation, and only 
the couple of e-folds around the end of inflation, reheating and the transition to 
radiation domination requires numerical treatment. Post-inflationary phenomena 
typically only span a few e-folds. 

Secondly, since we are doing quantum physics, the theory has to be renor- 
malized. In particular, the energy density which enters in the semi-classical 
Friedmann equation (see below) needs appropriate counterterms. Fortunately, 
features of the 2PI formalism include that it is renormalizable at any level of dia- 
gram truncation [211 [221 [231 [2H [25] , and that there is a similarly truncated energy 
density which is conserved. Hence by introducing (scale-factor dependent) coun- 
terterms for the energy density, mass and couplings, we can in principle cancel 
all divergences, and construct a well-behaved Friedmann equation. 
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We study a self-interacting real scalar, and go to NLO ((9(A^)) in a 2PI cou- 
pling expansion. By showing how to apply the procedure in practice, we expect 
it will be clear how to generalize to more complicated systems. In the conclusions 
we point out some issues, applications and ways of refining the approach. 



1.1 Setup 

We are concerned with a single scalar field with ip^ interaction. The action is 



S 
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(1) 



written in terms of co-moving spatial coordinates x. a{t) is the scale factor, and 
we assume a(0) = 1. Correspondingly, we will consider co-moving and physical 
momenta, denoted k and k = k/a(t), respectively. 

The evolution of the scale factor is in turn given by the Friedmann equation 
in terms of the Hubble rate 



H\t) 



(T°°(t)), 



3M2 



H{t) 



d{t) 



(2) 



Here we equate a classical quantity on the left-hand side to a quantum expectation 
value on the right-hand side. This only makes sense when the energy density is 
appropriately renormalized, an issue we will return to below. 

The system can be recast in comoving (conformal) time r], with dt = a{ri)dri 
and we can rescale the fielcil = (f){x) / a{ri) ^ in which case the action becomes 



S = I drjd^^ ^-{dr,4> - n4>) 



(3) 



We have introduced a new "comoving Hubble rate", Ti = a' /a = aH. In terms of 
the canonical momentum tt = d„(j) — Tid), the corresponding Friedmann equation 
isP 



/\2 



(4) 



In passing, it is useful to recall the classical equation of motion 



- dl 



a" A 

h aW + 

a 6 



</)(x) 



(5) 



^We use a to denote dta(t). a' will denote d^a(ri), with rj conformal time. 

^The equations of motion will be solved in terms of the "conformal" field 0(?7); but results 
converted back to the "physical" field ip{t). 

•^The right-hand side is not the Hamiltonian corresponding to the action but the Hamil- 
tonian of ([T]), written in terms of the rescaled fields. 
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The classical approximation amounts to generating a set of random initial condi- 
tions, solve for the evolution using ([5]) together with the Friedmann equation and 
then to average observables over initial conditions. In addition to using approxi- 
mate dynamics, also the classical averaging procedure is different from (jlj) in that 
the Hubble rate, and hence a{t), is derived from each individual initial condition 
rather than the average energy density. This could be resolved by simulating all 
initial conditions simultaneously using a common a{t) determined through the 
ensemble averaged energy density. Still, it would be a classical average rather 
than a quantum one. 

2 The 2PI formalism 

2.1 Equations of motion at 0{X^) 

We will not review the 2PI formalism in detail here, but refer to literature on the 
subject (for instance [26] and references therein). Applying the 2PI formalism to 
the action ([T]), while treating the scale factor a(t) as an external field, we find 
the equations of motion for the homogeneous mean field (f){ri) = {(pix)) and the 
propagatoiEl 

{T(j){x)(f){y)) - (t){ri)(t){i) = F{r], r]', x - y) - -p{r], r]', x - y) signc(r/ - r]'), (6) 

where F and p are real. 
The general form is 

[dl + e + M^'n)] F{r],r]',k) 

[d'^ + e + M'ir^)] p(r/,r/',k) 

These are known as the Kadanoff-Baym [27] or real-time Schwinger-Dyson equa- 
tions. They include a set of self-energies M^, M^, S,^, Si?, Sp, to be calculated and 
integrated up with the correlators for all past time. This makes 2PI simulations 
numerically challenging, although one should remember that 2PI simulations do 
not require statistical averaging, since the variables F, p, are already the full 
correlators. 

*The fields are defined along the Keldysli contour C. 



di j rf3xS<^(r/,V,x)0(r/'), (7) 



- / dr/%(r/,r/",k)F(r/",r/',k) 



cir/"S^(r/,r/",k)p(r/",V,k), (8) 



- / rfr/%(r/,r/",k)p(V',V,k). (9) 
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Throughout, we use Fourier transforms in comoving coordinates and assume 
homogeneity, so thatS 



x)0(r/', y)) = / (0k(r/)0-k(V))e^'^(''-^^ (10) 



k 



In the 2PI formahsm, the self-energies are determined from truncations of a 
(2PI-)diagram expansion. At NLO in a couphng expansion [2l|8], 

Mliv) = -^ + a\v)^' + ^{^\v)+mV,V,^ = 0)), (11) 

n"(n) A 
M'iv) = -^ + a'{r^)m' + -{^'{v)+F{v,V,^=0)), (12) 

}? - ( 1 \ - 

Si.(r/, V, x) = -—i){r]) \ F'^{r], v', x) - -p^(r/, v', x) j 0(r/') 

-^Fiv, V, x) (F'iv, V', x) - Ip\v, v', x)) , (13) 
Sp(r/,r/',x) = -y0(r/)2F(r/,r/',x)p(r/,r/',x)0(r/') 

V, x) hF^{r], 7]', x) - -p2(r/, r/', x) j , (14) 



A^ / 1 \ 

S</,(r/, V, x) = V, x) f 3F2(r/, 7]', x) - -p2(r/, r/', x) j 



(15) 



The only difference to the Minkowski-space equations is in the a(?7)-dependence 
of the effective masses. This is the standard free-field dependence, also appearing 
in the classical equation of motion Note however that the time coordinate is 
the conformal one and that we still consider the rescaled field (p. Physical time 
is therefore 

rv 

t -to= / a{T])dr], (16) 

and the physical correlators are 

(^(x)) = (^(x) = ^0(x), (17) 
a{r]) 

(TcpixMy)) = -i-^(T0(x)0(y)), (18) 

{TdMx)dMy)) = ^(TH-H)0(x)][(a,-H)0(2/)]). (19) 



^We use a unifying notation for momentum integrals in the continuum J^^ = J (2Tr}'i ^^'^ '-'^ 
the lattice J^ = y J2k- 
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For m = and with non-minimally "conformal" couphng to curvature (an addi- 
tional term \^'^R in the action), no trace remains of the expansion in the equation 
of motion, and physical processes proceed exactly as in Minkowski space, except 
that the time coordinate is "stretched" through ( |T6l) . We use minimal coupling 
to gravity and non-zero mass, leading to deviations from conformal behaviour. 



3 Renormalization 

Renormalization in an expanding space-time requires additional counterterms 
compared to the Minkowski operators involving the metric (in this case, 

the scale factor a{t)) emerge in the effective action. We will adopt the approach 
of [28], where the energy density is renormalized by subtracting a contribution 
corresponding to the adiabatic vacuum solution in the background defined by 
a{t). This vacuum solution can be solved for order by ordecl in derivatives of 
a{t) using a WKB-type ansatz: at leading order (d°) all quartic divergences are 
canceled, at the next order (a^) also quadratic divergences and finally at order 
the logarithmic divergences. 

In Minkowski space, a full-fledged 2PI renormalization procedure has been 
developed [211 ESI [23], [211 [25] , which allows a proper continuum limit to be taken. 
Since we are not intending to go to the continuum, we here take a somewhat 
simpler approach, and renormalize the equations of motion and the energy density 
in the LO (Hartree) approximation by a mass and energy counterterm only [8]. 
The mass counterterm is calculated using the same WKB vacuum solution as for 
the energy density. From a numerical viewpoint, logarithmic divergences are very 
small indeed, and so we restrict ourselves to renormalizing the energy density to 
this level of precision, i.e. expanding the WKB solution only to second order in 
d. 

3.1 Counterterms 

The energy density is, in the 2P1-L0 approximation {F is the correlator for the 
original 93 field) 

1 „ „ . 1 



2^2(i) <9x9x'F(t, t, X - x')x=x' + + t^^ = 0)) + 



24 



(^ip\t) + Q<f\t)F{t, t, X = 0) + 3F\t, t, X = 0)) - 5T°° .(20) 



^For brevity, we write ct" to mean all combinations of n time derivatives. For n = 2 we have 
terms and and a, etc. Similarly 0{H'^) includes 0{H). 
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We have introduced a bare mass mh{t) and an energy counterterm 6T^^(t) to be 
determined. These are time dependent, but only through a{t). The energy can 
be written in a more suggestive form as 

^ i (^dA'F{t, t', k),=,, + + + t, X = 0) + ^^2(0)^ t, k)^ 

t, X = 0) - ^F2(t, t, X = 0) - 5T00(t)] . (21) 
where ifQ = {ip{t = 0)). The LO equations of motion read 

m = 0, (22) 

[d^ + 3H{t)dt + Col{t)] F{t, t, k) = 0, (23) 

with 

ulit) = k'/a^t) +ml + ^F{t, t, X = 0) + ^^\t). (24) 

We will renormalize in a vacuum to be determined below. Assuming we have 
such a vacuum, let us define 

ml = m^- ^SF,^{t, t, X = 0), (25) 

where is a WKB approximation to the exact, time-dependent, vacuum 

solution Fyac{t, t, X = 0). 
Similarly, let us define 

5T°° = 6TZ + 6T^l (26) 
STZ = ljdA'SFUt,t',k),=t' + l(^-^^ + 

5T°? = -^^^^F.,e(t,t,x = 0)-^5Fl(t,t,x = 0). (28) 

We will be interested in evolving various, non-vacuum, initial conditions in 
real time, giving rise to some evolution of a{t). In principle, we could then also 
solve the vacuum equation numerically, simultaneously in the background of that 
same a(t), and self-consistently use the resulting F^ac for the counterterms of the 
non- vacuum simulation. 

Because the vacuum solution depends only on a{t) and its time derivatives, 
and because the renormalised LO equations of motion in vacuum are similar to 
the free ones (insert fl25l) into fl2^ intc l23l) ) . we will instead solve a WKB equation 
for the field modes. In this way, the counterterms will simply be functions of 
a{t), to be calculated at each time t. By including subsequent orders in WKB, 
we will be able to get a better and better approximation to -Fvac- 
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3.2 Choosing the vacuum 

Following [28], we have for the original, un-rescaled field ip a WKB ansatz 

¥^k(t)=ak/k(t) + aUW, m= nr4^^^ e-^^'''-^''^'^r (29) 

^J2a-\t)il]^{t) 

which is meant to satisfy the equation of motion for the free field modes </?k(0 

[dl + ?,H{t)dt + ulit)] Mt) = 0, (30) 

with 

Col{t) = e/a\t) + M\ = m^ + (31) 

As mentioned, we include the possibility of having a non-zero initial mean field 
<^(0), in which case we include the initial value in the mass. As is argued in the 
Appendix, divergences associated with a time-dependent mass (and in this case, 
a varying mean field) are logarithmic, and so beyond the level of approximation 
aimed at here. 

fik has to satisfy 

At leading order in WKB, fi^ = <jJk- By plugging this back into (!32|) we get the 
next order, 

a"/a ( lfaMV\ H'^ ( I f aMV 5 f aM 



atUkj \ 2 yacjk/ / (a^^kj \2 vatUk/ 4 yacuk 

where we have defined Ti = aH = a' /a. f^k defines our approximation to the 
(infinite order in WKB) vacuum. We define 

(5Fva.(t,t',x = 0) = Jj^l{t)ipk{t')). (34) 

Using 

^it)=(^-^n^-^H-^-^^Ut), (35) 

we have to this order and second order in H 

1 f ^ I ^ H? I ^ f 1 f m 



^ _y__a..ll + ^ll+^^j .-^^1 1|, (36) 



We also have 



5F.ac(t,t,x = 0)= / (37) 

determining the bare mass (l25l) and ST-^^ (l28l) . 

With m^, (JTf^gjj and 5TJ°^ as defined above, the renormahzed energy density 
has divergences oc C{t) In A in Minkowski space, where C{t) is a function of mass 
dimension four of the effective time-dependent mass and its time derivatives. In 
FRW space, C{t) can also depend on (time- derivatives of) H. In particular, 
in this approach the finite parts of the counterterms are chosen to cancel the 
1-loop corrections to 0{H'^), and so effectively amount to the renormalization 
conditions, in terms of some momentum cut-off A, 

(T°°)rcn = + C(C(0) In A), (initially in vacuum, LO), (38) 
m^^^ = + C»(C(0) In A), (initially in vacuum, LO). (39) 

Using the WKB solution with 0^ as an initial state, {T^^)^.^^ is identically 
zero for zero mean field, or initially equal to the "tree- level" energy density for a 
non-zero mean field 

(T°°(0))ren = [^(5*^)^(0) + ^mV(0) + ^^^(0)]. (40) 

In the simulations carried out here, H/uj\^ ^ 1/200, making the 0{H^) corrections 
very small indeed. Specifically including a non-zero cosmological constant is 
straightforward, but we will not do so here. 

Renormalization at NLO in 2PI in Minkowski space involves real-time solu- 
tion of separate auxiliary equations and counterterms for each included diagram. 
This is because the 2PI diagrams resum self-insertions to all orders, and so the 
structure of divergences becomes more involved. Generalization to expanding 
backgrounds becomes even more complicated and is beyond the scope of this pa- 
per. We will simply note that although formally quadratically divergent, vacuum 
corrections to the mass from the Sunset diagram are numerically more than an 
order of magnitude smaller than the LO contribution, also at the largest coupling 
used below, A = 6 [S]. At larger coupling this may no longer be the case, but 
then the coupling expansion should presumably be discarded altogether. Similar 
arguments apply to the vacuum contribution to the energy density. 

This concludes our treatment of renormalization. Obvious refinements are 
possible, in particular if one is interested in subtle issues like particle creation 
from the vacuum, thermal (or warm) infiation. But for our present purposes, 
this will suffice. 
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3.3 Initial conditions 

We are interested in two types of initial conditions. The initial correlators are 
chosen gaussian, of the form 



= '^in' . (dt^ldfP^.) = « + 1/2) D,in)Q^, (41) 
{[dtvlvi]) = h {{dt^l^k])=D2{n), (42) 
in terms of a dispersion relation 

Q"^ = Ok (M ^ Mgap) , (43) 
and a particle number in (approximate) equilibrium 



-1 



K= , Finite T, < = Jk^ + , (44) 



gap' 

and out of equilibrium 

= c, |k| < c, = 0, |k| > c, Step. (45) 

T, c are free to be chosen. T = is the vacuum. The expansion is encoded in f^k 
and the correction factors 

^ 1 . 7^ f 1 4%V . i (%VV -^(^A (%V 1 .(«) 



(<)2 \^ V < / 4 V < / y < \ ■ ■ 2 V cuj- 

Mgg^p is found by solving the 1-loop gap equation for the relevant initial condition, 

Ml, = mU\ {m + Fo) , Fo = J^ ^^i^i^. (47) 

The "Step" represents some generic out-of-equilibrium initial state, which is vac- 
uum for large k. More elaborate choices are of course possible. "Finite T" is 
strictly speaking only equilibrium in the Hartree approximation at if = and at 
the initial time. 



3.4 Observables 

We monitor the global quantities physical time f|T6l) . scale factor a{t), the Hubble 
rate, the mean field 0(t)/Mpi, renormalized local correlator 

F,,^{t, t, X = 0) = F{t, t, x = 0) - (5F,ac(t, t, X = 0), (48) 

and renormalized energy density (T°'^(t))i.cn- 
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3.4.1 Particle numbers 

For consistency witli tlie renormalization and initialization prescriptions above, 
we should define particle number relative to the vacuum defined by Cl^ (l33i) . to 
the same order in H. We use 

By inserting 0^, we find 

,3(,),..„) ^^_l,J^Ji)\lhJf\\^o . (50) 

As expected this is zero in Minkowski space. The particle number relative to our 
vacuum is therefore 

Out of equilibrium, it is advantageous to use instead the "self-consistent" defini- 
tions [291 EDI E] 

a'(tK(t) = a3(t)v/(0k(t)0k(t))(^k(t)<^k(t)) -\- a'(t)nr(t), (52) 

uf = V{Mt)Mt))/{Mt)Mt))- (53) 

Way out of equilibrium, lj^ will look wild, but close enough to equilibrium the 
two definitions agree in Minkowski space. With expansion it is easy to see that 
in our vacuum nk(t) = nk{t) + O (a^), and that 

,54) 

In practice, the two definitions of n\^{t) were seen to agree. 

In the context of thermalisation, we quote the effective temperature and chem- 
ical potential of a mode, based on assuming a Bose-Einstein thermal distribution 
asymptoticalljlll, 

, .off ,,off\ / 1 \ , ,cff i.cff 

^k ~ /^ch \ 1 \ , 1 / , 1 \ ^k ~ /^ch 



When and are independent of k, we say that the system has equilibrated 
kinetically. In Minkowski space the timescale for this to happen for A = 6 (as 
used throughout most of this paper) is ~ lOOOm"^ |8]. Chemical equilibration 
is when this common is zero. Again in Minkowksi space, this is roughly an 
order of magnitude slower than kinetic equilibration. 



''We note that a^(t)rik rather than rik(i) itself is the quantity that is expected to equihbrate 
to a Bose-Einstein distribution, as a result of defining Fourier transforms in terms of co-moving 
k. 
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4 Example applications 



The setup introduced above can now be applied to systems and phenomena of 
interest. We use lattices of 32^ points, time-steps 5t = 0.05, lattice spacing 6x = 
1. The masses are m/Mp\ = (massless) m/Mpi = 0.00025 — 0.001 (massive), 
and for the interacting case A = 6 to study thermalisation, A = 0.1 for the case 
of preheating. Thermal initial conditions have T/Mpi = 0.0005 — 0.002 and out- 
of-equilibrium Step initial conditions have c/m = 5. We use Mpi to set the scale 
of gravity and hence the expansion rate. Large Mpi means slower expansion. 
When a{t)m ^ 1, one should start worrying whether the lattice is too coarse. 
When a{t)m > 2, results can probably no longer be trusted. In the massive case, 
a(0)m = 0.2. 

From the point of view of a realistic cosmology, at least after inflation our 
choice of relative Mpi is much too small, and our coupling too large. In m'^ip'^ 
inflation m/Mpi < 10^^ and in A<y9^ inflation A ~ 10"^^, to be consistent with 
measurements of the Cosmic Microwave Background. Typical reheating temper- 
atures are or the order of 10~^Mpi. Therefore the parameters here over-emphasize 
expansion effects compared to most post-inflationary cosmological phenomena. 
Inflation itself is an exception, since the relatively much larger couplings used 
here (in the preheating application, 0.1 instead of 10~^^) may make larger inter- 
action rates compensate for the larger expansion. 

It may also be worth pointing out that as H < 10~^ (see below), all the modes 
under consideration are sub-horizon. 



since we only allow a(t) < 10. 

4.1 Free field at finite temperature 

The simplest case is a free field initially in the vacuum. But since we have 
effectively chosen zero cosmological constant as our renormalization condition, 
the Universe is just static. Still, it is a useful test of the numerics and the 
renormalization. 

A slightly more interesting test case is to initialize a free field in a thermal 
state, A = 0, T/Mpi = 0.04, m/Mpi = 0.02; 0, shown in Figure [H The effect of 
the expansion is to redshift the momenta kphys(t) = k = k/a(t), but the particle 
number as a function of co-moving momentum a'^nk(t) is constant. If the field is 
massless, we have 



with in our case f{x) = (exp(x) — 1)~^. If we identify the constant T as the 
temperature, then /(|k|/T) is constant if T{t) = T(0)/a(t). Hence for the mass- 
less case, T^^ (l55l) should be independent of k and decrease as 1/a. These are 



^ 0.2/a(t) > H, 



(56) 



aV(t) = /(^k/T) = /(|k|/T), 



(57) 
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k ^ /M , 

phys pi 

Figure 1: The effective temperatures T^^ for all modes /c in a free simulation 
starting from equilibrium. In the massles case (black), the modes stay thermal 
with a temperature dropping as 1/a. In the massive case (red) this is only true 
in the UV. The range of physical momenta redshift towards the IR. Overlaid are 
(dashed) lines denoting T{t = 0)/a{t) in the two cases. 



the black lines of Figure [H In the massive case, we no longer have (1571) . since 
cDk = \/rn?~+~\k^Ja? . For |k| ^ m the argument still holds, but for small |k| 
we have deviation from |k|-independence. This is shown by red lines. We may 
choose to interpret the deviation as a chemical potential, and it is easy to seqj 
that for |k| ^ m, T oc l/a^(t) and fil^/m = 1 — l/a^(t). We stress that at this 
point the chemical potential is only meant to express a relative over- abundance 
of particles in the low momentum modes. 

As a consequence of redshift only, we therefore expect that a non-interacting 
massive field will have an effective chemical potential going asymptotically to 
fi'^^/m = 1. When including interactions, these will drive the chemical potential 
towards zero, but only if they are fast enough on the time-scale of the expansion. 

We can trust the simulation as long as the high |k| modes stay in the vacuum, 
with no significant thermal population. This is automatic in the free field case, 
but when including interactions, high |k| modes will be excited. 
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Figure 2: The scale factor (upper left), mean field (lower left), renormalized 
correlator and energy density (right) in the presence of an oscillating mean field. 
Both correlator and energy density are well renormalized, and the energy density 
behaves as matter, oc a~^. The scale factor evolves as matter domination a = 
(1 + t/to)^'^^ (overlaid), and the mean field amplitude decreases as 1/t (overlaid). 



4.2 Free mean field and vacuum modes 

Another interesting check of the numerics is to initialize the mean field away 
from zero and let it oscillate freely, without interacting with the modes. These 
are initialized in the vacuum. It is easy to see that an oscillating, homogeneous 
mean field behaves as matter, with zero pressure, {T^'^)j.cn oc and a{t) oc t'^^^. 
Also, the mean field is expected to have the form 

^(^) = TT?7^'°'^"'^^- ^^^^ 

This supposes that the vacuum has been correctly renormalized and does not 
contribute to the energy density. In principle, fast expansion could induce particle 
production (an 0{H^) effect), but this will be very small indeed. 

In Figure [2] we show the scale factor (upper left), the mean field (bottom left) 
and the renormalized energy density (bottom right). They all scale with a{t) as 



^Use that Wk ~ m + k'^/{2a'^m). Then for pk - /Uch(0]/r(t) to be constant, ^°^(t) 
m(l - l/a'^{t)). 
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Figure 3: Left: The mode spectrum in time, starting from a thermal initial 
condition. The observable is chosen so that it is a straight line in kinetic equi- 
librium fl55l) . and effective temperature and chemical potentials can be found 
from a fit (overlaid). Right: Effective temperature and chemical potential for 
a fast (r"^'VMpi = 0.008), slow (T'^^VMpi = 0.002) and very slow expan- 
sion (T™'*/Mpi = 0.0005). The temperature drops as matter domination, and 
the chemical potential becomes larger for faster expansion, when interactions 
have a harder time keeping up. Dashed lines represent the free field behaviour 
(Ji'ch/m = 1 - l/a?{t). 



expected. In the top right frame, we show the renormalized energy density after 
the mean field contribution is subtracted, and the renormalized equal time cor- 
relator. Both are renormalized to 10~^ (here in units where the un-renormalized 
quantity is 0{1)). 

4.3 Interacting field at finite temperature 

Interactions alter the dynamics and the thermal state is only approximately Bose- 
Einstein, with an additional thermal mass component. A fair approximation to 
an initial thermal state is to solve for the thermal mass at LO, and use this to gen- 
erate a Bose-Einstein distribution (jH]). Because we evolve with NLO equations 
of motion, the initial condition is not exactly the equilibrium one. Also, since we 
include expansion, the system can anyway only be in approximate equilibrium. 
We use A = 6, Ti„it/Mpi = 0.0005, 0.002, 0.008, T;^it/m = 2. 

Figure [3] (left) shows ln(l/nk + 1) (l55ll in time. We see that interactions 
partially compensate for the chemical potential /ich/"^ = 1 — l/'2^(^)) although 
kinetic equilibration is not complete. We may tentatively extract effective tem- 
peratures and chemical potentials from fits to the spectrum, shown in the right- 
hand panel. In all cases, temperature drops roughly as matter domination. As 
for the non-interacting case, the chemical potential rises and goes asymptotically 
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to a finite value. Clearly, interactions are not fast enougli to upliold or restore 
cliemical equilibrium. In the fast expansion case, although the rise is slower than 
m(l — l/a^(t)) (the initial mass is not the zero-temperature one), asymptotically 
it goes to m, suggesting that interactions are ineffective. In the slow case, the 
asymptotic value is about 3/4 as large, as a result of interactions. In the very 
slow slight decrease can be seen at late times. Still, in all relic 

abundance of particles freezes in. 

A common criterion for staying thermal is for some interaction rate F to 
dominate the Hubble rate. Such a rate could be provided by the Sunset diagram 
damping rate [HI]. At large temperature T ^ 



The Hubble rate is given in terms of the energy density 



r(2^,A)^^^. (59) 



= T^ (60) 



and so with our parameters 



^ ^ 0.09^ » 1, (61) 

suggesting that the system should be able to stay thermal. Even more so as 
temperature drops, although when m ~ T flM]) no longer applies. The "scattering 
time scale" is mtgcat. — — 17. Chemical equilibration is known to require 
of order SOOtgcat. at this coupling in Minkowski space [8]. Therefore T/H ^ 1 
is not sufficient for chemical equilibration, although as we will see below, kinetic 
equilibration does take place. 

The 2PI-resummed Sunset diagram includes 2-to-2 scattering as well as off- 
shell l-to-3 scattering. The damping rate fl59l) is dominated by the former, which 
leads to kinetic but not chemical equilibration, as particle number is conserved. 
The latter off-shell process does change particle number, but is a higher order 
effect included in the diagram through the 2PI resummation of many perturba- 
tive diagrams. It is therefore no surprise that chemical equilibration happens 
on a much longer timescale, being (naively) suppressed by additional powers of 
the coupling. This highlights the possible shortcomings of criteria like (16T]) to 
establish thermalization. 



^We note that our system is probably not in this asymptotic regime, and not at small 
coupling, and so a more complicated expression may be required for detailed estimates. For an 
order of magnitude estimate, however, this will do (see also [8]). 
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Figure 4: The spectrum a^nk vs. k when starting from a Step initial condition. 
Again, left is for T"^'7Mpi = 0.002, right for T"^'7Mpi = 0.008. Kinetic equili- 
bration happens on timescales of order mt = 200 in both cases. The resulting 
temperature is clearly smaller for fast expansion. 

4.4 Thermalisation with expansion 

In order to study kinetic equilibration, and its dependence on expansion rate, we 
start the system way out of equilibrium using the Step initial condition. Figure H] 
shows the evolution of the spectrum at early times, for a slow expansion (m/Mpi = 
0.001, left) and a fast one (m/Mpi = 0.004, right). The Step is smeared out into 
a smooth spectrum on a timescale of mt ^ 200, while the modes are redshifted 
towards the infrared. In this particular case, redshift may help equilibration, as 
evolution from the Step to equilibrium involves transfer of power to smaller |k| 
modes. For even faster expansion, this redshift may "overtake" the equilibrium 
distribution, and intermediate- |k| modes will have to be re-populated through 
scattering. 

Once the spectrum is kinetically equilibrated, evolution is much slower, and 
we can again quantify the spectrum by the effective temperature and chemical 
potential. This is shown in Figure O with the on the left and the fi^^^^ 
on the right. Temperature drops as with a between 0.4 and 0.8, while the 
chemical potential again asymptotes to a finite value. Faster expansion (smaller 
Mpi) again results in a larger asymptotic value. 

We are therefore in a regime where interactions are strong enough to equi- 
librate an initial out-of-equilibrium condition into a Bose-Einstein-like thermal 
state, but not strong enough to get rid of the chemical potential. As in the sec- 
tion above, we conclude that the T/H ^ 1 criterion applies to restoration and/or 
maintaining of kinetic equilibrium only. 



17 



1.1- 



....Mp=2000m 
....Mp=1000m 
....Mp=500m 
....Mp=250m 



1.5 



1 - ■ 



0.5- 



T/m 



,. ..M =2000 m 
pi 

....Mp = 1000m 

,...Mp=500m 

....M =250 m 
pi 



"■^0 1000 2000 3000 4000 5000 6000 ''o 1000 2000 3000 4000 5000 6000 

mt mt 



Figure 5: Effective chemical potentials (left) and temperatures (right) in time 
for different expansion rates. Only after some time are the spectra sufficiently 
straight to allow fits to determine T'^^ and nl^. The asymptotic chemical potential 
increases with increasing expansion rate, interactions become more inefficient and 
more particles freeze in. Final temperature decreases with expansion rate as a 
power law, ^-(O-^-o.s)^ 



4.5 Preheating after inflation 

At the end of cosmological inflation, the inflaton mean field leaves the slow- 
rolling stage and begins oscillating around a minimum of its potential. Through 
interaction with other fields and/or its own field modes, energy is transferred into 
particle excitations of these fields. In some cases, resonant particle creation can 
take place, known as preheating [32]. This usually lasts for a few periods of the 
inflaton oscillation, during which some fraction of the energy is transferred, after 
which normal perturbative decay and reheating transfers the rest. 

Let us assume for simplicity that the mean field oscillates harmonically as 
( I58p . including the self- interaction (which induces an-harmonic oscillations) only 
through an altered frequency mosc 

4>(t) = (t)o(t) cos(mosci)- (62) 

Then each mode function ipk{t) obeys approximately 

V3k + 3H{t)ip^ + —— + + -(f)l{t) cos^(mosct) U'k = 0. (63) 
Wit) 2 J 

Again approximately, a mode k will be in resonance if 

= ^T-T — H — + 2 ~ ^ ' i^^> 



osc 



with / an integer, with the strongest resonance for = 1. In the resonance, par- 
ticle number grows exponentially. Because a{t) and the amplitude (f)o{t) decrease 
in time, in an expanding background, modes will move in and out of resonance. 
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Figure 6: An oscillating mean field coupled to its modes results in reheating and 
possibly resonant preheating. Shown is the mean field (upper) and the energy 
components (lower) for three cases: Large field, slow expansion (left); large field, 
fast expansion (middle), small field, slow expansion (right). For the first case, 
around mt = 100 the mean field dumps all its energy into particles through 
resonant preheating. 



Figure El upper panels, shows the mean field for three different cases, nor- 
mahzed to their initial amplitude: m/Mpi = 0.001, (^(0)/Mpi = 0.01 (black); 
m/Mpi = 0.004, <^(0)/Mpi = 0.04 (red); and m/Mpi = 0.004, <^(0)/Mpi = 0.01 
(red). Since energy density is dominated by the mean field, the expansion rates 
are the same initially in the first and second case, faster in the third. But Mpi 
is different in lattice units and in units of m in the first case. Around time 
mt = 100 the mean field oscillation suddenly collapses for the slow expansion, 
small m/Mpi, case only. In the lower panels, we see that at the same time, energy 
is transferred from the mean field component {E{(j))) to the mode component 
{E{G)) (normalised to the total energy). We conclude that we have preheating 
only in this one case. 

Indeed, in figured (left), we show particle numbers of the 5 lowest momentum 
modes; these grow exponentially, with the first (zero) mode leading the way early 
on reaching rik — 10000. But also in the fast expanding case (red), the zero 
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Figure 7: Left: Particle numbers of the 5 lowest |k| modes. Only in the large 
field cases is there a resonance, and only in the slowly expanding cases (black) 
is the resonance long enough to get large particle numbers. With fast expansion 
(red) only the zero mode grows briefly. Right: for the zero modes flMl) . with 
the replacement 0o ~^ 4'{t)- Hence proper is the envelope of the oscillating 
curves. 



mode grows. In that case the resonance ends rather early on, particle number 
is only order 100, and the non-zero modes are not excited. In the right panel, 
y4k for the zero mode is the envelope of the oscillating curves, shown for the 
three cases. The overall normalisation will depend on mosc, and apparently the 
resonance band is somewhere above Ak = 1. The qualitative picture is fairly 
clear. The blue curve is below the resonance, and no preheating occurs. The 
(amplitude of the) black curve goes through the resonance band long enough for 
a significant amplification of the zero mode. This eventually excites the non-zero, 
non-resonating modes to complete reheating altogether. For faster expansion 
(red), the mean field amplitude decays very fast, and apparently the zero mode 
is only briefly in resonance, but does not grow sufficiently to trigger substantial 
preheating. 

Clearly a more detailed study is required, for finer lattices where not only the 
zero mode resonates, and taking into account that the mean field oscillation is 
anharmonic. Also, in most realistic models of inflation, the mean field starts out 
0{Mpi). Finally, a realistic (self-) coupling would be many orders of magnitude 
smaller than the one employed here, A = 0.1. Such a study is beyond the scope 
of the present work. Still, resonant preheating is within the range of application 
of the present formalism, and outcomes depend strongly on the expansion rate. 

Parametric resonance was studied using 2PI methods in [33], in Minkowski 
space. Particle numbers grow exponentially to nk — 1/A, and so in principle a 
2PI coupling expansion like the one used here may not be applicable. Indeed, for 
small coupling, numerical instabilities are experienced. In [33] a 1/A^ expansion 
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in the number of fields is used, making the dynamics more stable. Using the 
procedure outlined here, this can be generalized to FRW spaces. 



5 Conclusion 

We have seen that extending the 2PI formalism of out-of-equilibrium quantum 
fields to expanding backgrounds amounts to introducing a time dependent mass 
a^m? — a" /a in the conformal time, rescaled field, equations of motion (10- 
19). Having solved these equations on the lattice, observables are translated back 
to physical fields and physical time. The scale factor is derived from the semi- 
classical Friedmann equation involving the renormalized energy (T°°)ren- We have 
here opted for an approximate renormalization strategy, where counterterms for 
the mass and the energy density are calculated in the LO/Hartree approximation 
in terms of a particular vacuum, the adiabatic free-field solution to second order, 
both in WKB and in a. We argued that going beyond this order is possible, 
although only really necessary when aiming at taking the continuum limit or 
using very large expansion rates. A fully 2PI renormalization beyond LO is much 
harder. 

On a finite comoving lattice, modes are redshifted towards the IR. Therefore 
the number of e-folds of expansion available is limited; a simulation can only be 
trusted as long as there is a range of UV modes that stay in the vacuum. Other- 
wise cut-off effects will influence the physics and presumably the renormalization. 
One way of quantifying this is for the time-dependent mass am to stay less than 
unity. In the "conformal" case m = another mass scale (temperature, initial 
mean field) will play a similar role. 

At the end of the day this is a practical question of computer capacity. 2PI 
simulation are memory intensive in that the memory of past time-steps must be 
saved to generate the self-energy kernels (right hand sides of (10-12)). On the 
other hand, no statistical averaging is necessary as the solutions to the equation 
are the full correlators. The total simulations performed here amount to about 
5000 CPU hours. 

Possible applications range over all the topics already studied in Minkowski 

3 El El El E], which enters in early Universe physics as 
reheating and preheating, both the resonant variety 
[33] and tachyonic preheating [Ml [35] . In the latter cases, it may be prudent to 
use a diagram expansion in 1/A^ rather than A, to be sure corrections are under 
control. 

In the present paper we have made test-runs of many of these cases and 
pointed out the main effects of the cosmological expansion. All of these applica- 
tions deserve further scrutiny, also of combined effects of smaller /larger couplings, 
expansion rate, temperature and masses. A study of the inflationary regime was 



space: Thermalisation [21 
well as heavy ion collisions 



■'in heavy ion collisions the expansion is somewhat different from FRW. 
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not attempted, mainly because it would require very large lattices (lots of expan- 
sion) and/or m = which is a very special case. 

One interesting result is the difference between the kinetic and chemical equi- 
libration timescales, and the possibility of having one but not the other happen. 
Clearly, one must be careful when assuming instantaneous thermalisation, as is 
sometimes done when considering reheating, preheating and phase transitions in 
the early Universe. Using a criterion like T/H ^ 1 to ensure instant thermalisa- 
tion presumes careful consideration of which F is the relevant one. 

We were also able to confirm that the amount of resonant preheating depends 
sensitively on the rate of expansion. This is because field modes are redhifted 
in and out of resonance bands. This shortens the resonance time, making even 
exponential growth much less effective. 

The quantum 2PI equations have a classical counterpart, i.e. 2PI-truncated 
equations for classical correlators, reproducing classical dynamics [3S]- The rela- 
tion between the two amounts to neglecting terms like compared to terms like 

in the self-energies (14-16), and ignoring renormalization. Although the full 
classical approximation is in principle exact, it relies on statistical averaging over 
initial conditions. Classical 2PI has no statistical errors, but diagram expansion 
truncation introduces a different type of approximation. As such, it contributes 
an alternative way of doing classical simulations. 

Extension to more complicated models than a single self-interacting scalar 
is straightforward. Possible applications include reheating and preheating with 
multiple fields (including fermions [5]), departure from and return to equilibrium 
for systems with heavy particles decaying into light ones. In the context of multi- 
field preheating this may allow to calculate non-gaussian signatures in the CMB 
[23 [M]. 

In conclusions, we believe that given sufficient numerical capacity, and observ- 
ing certain simple rules, the 2PI formalism for out-of-equilibrium fields provides 
a convenient tool for quantitative calculations in cosmology. 
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A Structure of divergences 

Here, we present a detailed exposition of the cancellation of divergences in the 
renormalization procedure presented in the main text. 
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We begin in Minkowski space, at 2PI-L0. At this order, the equation of motion 
for the propagator modes reads, 

[d^ + ul{t)]F{t,t,k)=0, (65) 

with 

uUt) =k' + ml + t, X = 0) + ^^\t). (66) 

Away from the vacuum, and for the moment neglecting time derivatives of the 
mass, writj"1 

Fit, t, X = 0) = F^ac + Fres(t) = / ""^^^.2^^^ (67) 



In our prescription, we choose 



F^ac = / , (68) 



and renormahze the mass by 



ml + = m^. (69) 



We will also need 



Then 



(70) 



Olit) = k^ + m^ + ^Fres(t) + ^<^'(t). (71) 



The energy density can be written 

(T°°(t)),e„ = lj^{dA'+^l{t))F{t,t\k),=, 



-^F(t, t, X = 0)^ + ^^\t) + A^^(t) + 5T°°. (72) 



We can rewrite this as 



yj^(nk + l/2)u^{t) - + - ^F^,, (F,es + ^'(t))(73) 



^^Note that although we have an interacting theory, at 2PI-L0 the field can still be written in 
terms of mode functions, and are therefore completely described in terms of a particle number 
and a dispersion relation with an effective time-dependent mass. 
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The first line has at most logarithmic divergences. The second line quartic, 
quadratic and logarithmic. Let us assume that the particle number decays 
faster than k~'^, in which case the divergencies result only from 

j^^ = F, + -^ {^\t) -^1 + F^t)) Fvac + O (logs) . (74) 

The middle term almost precisely cancels the last term in (1731) . leaving only a 
divergent constant. Therefore, by choosing 

ST'' = + ^^IF^^, - Fo, (75) 

all the remaining quartic and quadratic divergences are cancelled. 

Let us now consider what happens when including time derivatives of the 
mass. We can consider this in the WKB approximation, in a similar way to the 
main text, but with 

ul = e + M\t), M^=ml + ^ + Fvac + F,es(t)) . (76) 
The ansatz is 

Mt)=aMt) + aim, Mt) = -^=e-^I'''-^''^'r (77) 

V ^i'k(rj 



As above, we find that 

2 



4U 



cjk\ MM f 



with 



k / Minkowski k \ k / Minkowski 

We can compare this to the FRW case discussed above, where we have 



When we calculate the energy density, we have for instance 

2 



fik / , , 1 [ ^ 



(7J 



(79) 



-i = --(^-^]- (80) 
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The leading term in time derivatives (no derivatives) is the usual quartically 
divergent term. In the FRW case, the next order is of the form 



i 



or H 



(82) 



oc 



which is quadratically divergent. However, from the time derivative of the mass, 
we get terms of the form 



which are logarithmically divergent as well as finite terms, and in FRW we in 
addition get mixed terms like 



which are also logarithmically divergent. 

Therefore, since we choose to ignore logarithmic divergences in the energy, no 
counterterms proportional to the derivative of the mass are necessary. 
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